
clear all;
clc;
disp('forecast one-quarter ahead');
load ('optimization_g1000k_z(pd)_40qtrs_.mat','NUMBER_OF_YEAR','Quarter_LIST_LENGTH','RATE_LIST_LENGTH','CRITERIA','R','C','N','stdProbQ','defaultFreq','bin','avgRollingTM','naive1','naive2');

% INSAMPLE_YR = 2008; % !!!!
INSAMPLE_YR = 2007;
INSAMPLE_QR = 4;


inSampleLength = (INSAMPLE_YR - 1998)* 4 + INSAMPLE_QR;
outSampleLength = NUMBER_OF_YEAR * Quarter_LIST_LENGTH - inSampleLength;



% To calculate credit cycle index Z, real and forecast----------------------
% To minimize the distance to get estimate of w----------------------------

% Website tutoring upload data from xls to matlab:
% http://blinkdagger.com/matlab/matlab-using-xlsread-to-import-excel-data/

% [fcstInvPD] = xlsread('fittedReg_SampleEndTo_07Q4_09Q3_PD.xlsx','s51:s54'); %insample 44qtrs !!!
[fcstInvPD] = xlsread('fittedReg_SampleEndTo_07Q4_09Q3_PD.xlsx','s47:s54');%insample 40qtrs


P = stdProbQ;
Z_fcst = zeros (outSampleLength, 1);
Wfcst = zeros (outSampleLength, CRITERIA);
Pfcst = zeros(CRITERIA, outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
efcst = zeros(outSampleLength,CRITERIA,3); % 3 columns means model, avg and previous methodology
avgPDfcst = zeros(outSampleLength,4); % avg default prob for model, benchmarks and empirical
down = zeros (CRITERIA,outSampleLength,4);

lb = zeros(1,inSampleLength);
ub = ones(1,inSampleLength);
A = []; b =[];
Aeq =[];beq=[];
x0 = 0.2;


for i = 1: outSampleLength
    % To calculate credit cycle index Z, real and forecast----------------
    
    t = inSampleLength + i -1; % t is 44 for calc w for fcst 45th Q
    PD = defaultFreq(1:t);
    inversePD = norminv(PD);
    Z = - zscore(inversePD);
%     inversePD(t+1) = fcstInvPD(i); % new try for including fcst PD to calculate Zfcst
    Z_fcst(i) =  - ( fcstInvPD(i) - mean(inversePD) ) / std(inversePD);
    
    
    % To minimize the distance to get estimate of w------------------------
    % Using historical 44 qtrs to estimate w, so we can use it to fcst 45th
    
    X(:,:) = bin(t,:,:);
    
%     % L1
%     obj1=@(wfcst)difL1_MultiObj(wfcst,t,P,X,Z,R,C);
%     [Wfcst(i,1),fval] = fminimax(obj1,x0,A,b,Aeq,beq,lb,ub);
%     
%     
%     % L2
%     obj2=@(wfcst)difL2_MultiObj(wfcst,t,N,X,Z,R,C);
%     [Wfcst(i,2) ,fval] = fminimax(obj2,x0,A,b,Aeq,beq,lb,ub);
%     
%     
%     % WAD
%     obj3=@(wfcst)difWAD_MultiObj(wfcst,t,P,X,Z,R,C); % weight is est. P
%     [Wfcst(i,3),fval] = fminimax(obj3,x0,A,b,Aeq,beq,lb,ub);
%     
%     
%     % NAD
%     obj4=@(wfcst)difNAD_MultiObj(wfcst,t,P,X,Z,R,C);
%     [Wfcst(i,4),fval] = fminimax(obj4,x0,A,b,Aeq,beq,lb,ub);
%     
%     
    % SVD
    obj5=@(wfcst)difSVD_MultiObj(wfcst,t,P,X,Z,R,C);
    [Wfcst(i,5),fval] = fminimax(obj5,x0,A,b,Aeq,beq,lb,ub);
    
    
%     % NSD distance
%     obj6=@(wfcst)difNSD_MultiObj(wfcst,t,P,X,Z,R,C);
%     [Wfcst(i,6),fval] = fminimax(obj6,x0,A,b,Aeq,beq,lb,ub);
%     
%     
%     % BFS distance
%     obj7=@(wfcst)difBFS_MultiObj(wfcst,t,N,P,X,Z,R,C);
%     [Wfcst(i,7) ,fval] = fminimax(obj7,x0,A,b,Aeq,beq,lb,ub);
%         
%     
    % D3
    obj8=@(wfcst)difD3_MultiObj(wfcst,t,P,X,Z,R,C); 
    [Wfcst(i,8),fval] = fminimax(obj8,x0,A,b,Aeq,beq,lb,ub);
    
    
    % D1
    obj9=@(wfcst)difD1_MultiObj(wfcst,t,P,X,Z,R,C);
    [Wfcst(i,9),fval] = fminimax(obj9,x0,A,b,Aeq,beq,lb,ub);
    
    % D1^2
    obj10 = @(wfcst)difD1sqr_MultiObj(wfcst,t,P,X,Z,R,C); 
    [Wfcst(i,10),fval] = fminimax(obj10,x0,A,b,Aeq,beq,lb,ub);
    % subplot(5,2,8); fplot(obj8,[0 1],'*');
    % title('Mininization Object with D3 distance',...
    %   'FontWeight','bold');
    
    
    
    %To calculate out-of-sample fcst matrix--------------------------------
    for type = 1: 10
    Z(inSampleLength+i) = Z_fcst(i);
    Ptemp= calculatePHat_fixedW(Wfcst(i,type),t+1,X,Z,R,C); % PHat = f(X,What,Zhat) and t+1 to calc fcst 45th Q
    Pfcst(type,i,:,:) = Ptemp (t+1,:,:);
    clear Ptemp
    end
    
    
end


% Forecast accuracy--------------------------------------------------------

for i = 1: outSampleLength
    
    % error(Mean Absolute Error)
    
    t = inSampleLength + i -1;    
    x2(:,:) = avgRollingTM (t,:,:); % take average as naive1 benchmark
    x3(:,:) = stdProbQ (t,:,:);% take previous as naive2 benchmark
    y(:,:) = stdProbQ(t+1,:,:);
    num(:,:)= N(t+1,:,:);
        
    for type = 1: CRITERIA
        x1(:,:) = Pfcst(type,i,:,:);    
        efcst(i,type,:) = calculateMAE(type,x1,x2,x3,y,R,C);
    
    
      
    % PD comparison between model and benchmarks.
    % 1. Average PD along time 
    avgPDfcst (type,i,1) = mean(x1(:,C));  
    avgPDfcst (type,i,2) = mean(x2(:,C));
    avgPDfcst (type,i,3) = mean(x3(:,C));
    avgPDfcst (type,i,4) = mean(y(:,C));
    
%     % Below method cannot give u good result
%     avgPDfcst (type,i,1) = sum(x1(:,C).* sum(num,2))/sum(sum(num)) ;  
%     avgPDfcst (type,i,2) = sum(x2(:,C).* sum(num,2))/sum(sum(num)) ;
%     avgPDfcst (type,i,3) = sum(x3(:,C).* sum(num,2))/sum(sum(num)) ; 
%     avgPDfcst (type,i,4) = sum(y(:,C).* sum(num,2))/sum(sum(num)) ; 
        
    % Downgrade 
    
    for ii = 1: R
        col = ii+1;
        for jj =col: C
            down(type,i,1) = down(type,i,1)+ x1(ii,jj);
            down(type,i,2) = down(type,i,2)+ x2(ii,jj);
            down(type,i,3) = down(type,i,3)+ x3(ii,jj);
            down(type,i,4) = down(type,i,4)+ y(ii,jj);
        end
    end
    
    end
end


cellProb5_outsample = cellProbCompare(Pfcst,naive1, naive2, stdProbQ, R, C, outSampleLength,5,1);
cellProb10_outsample = cellProbCompare(Pfcst,naive1, naive2, stdProbQ, R, C, outSampleLength,10,1);

%Dif criteria Default Prob (10 = D1^2; 5 = SVD)
figure
for i = 1:R
    temp(:,:) = cellProb10_outsample(i,C,:,:);
    subplot(outSampleLength/2,2,i); plot (temp);
    title('In-sample D1^2 fitted Default Probability Comparision',...
        'FontWeight','bold')
    set(gca,'XTick',1:4:outSampleLength);
    labels = quaterlabels(INSAMPLE_YR+1, outSampleLength);
    labels = labels(1:4:outSampleLength);
    set(gca,'XTickLabel',labels);
    xlabel('Quarter');
end
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');

% Dif criteria diagonal Prob

clear temp
figure
for i = 1:R
    temp(:,:) = cellProb10_outsample(i,i,:,:);
    subplot(outSampleLength/2,2,i); plot (temp);
    title('In-sample D1^2 fitted Diagonal Probability Comparision',...
        'FontWeight','bold')
    set(gca,'XTick',1:4:outSampleLength);
    labels = quaterlabels(INSAMPLE_YR+1, outSampleLength);
    labels = labels(1:4:outSampleLength);
    set(gca,'XTickLabel',labels);
    xlabel('Quarter');
end
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');

figure
f(:,:) = avgPDfcst(10,:,:);
plot (f);
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');
set(gca,'XTick',1:1:outSampleLength);
set(gca, 'XTickLabel', {'08Q1','08Q2','08Q3','08Q4','09Q1','09Q2','09Q3','09Q4'});
title('average Rating Default Probability Forecast: 2008Q1 - 2009Q4');

%  % 2. PD along rating
%     pdQ (:,1) = x1(:,C); pdQ (:,2) = x2(:,C); pdQ (:,3) = x3(:,C); pdQ (:,4) = y(:,C); %each rate default in every quarter
%     
%     subplot (outSampleLength/2,2,i);plot (pdQ);
%     legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');
%     title('Default Probability Forecast for each Q: 2008Q1 - 2009Q4');
%     

% D1sqr metrics for model and benchmarks
%     for ii = 1: RATE_LIST_LENGTH-1
%         for jj = 1: RATE_LIST_LENGTH
%             s(i,1) = s(i,1) + (ii - jj)* sign(x1(ii,jj) - y(ii,jj))*(x1(ii,jj) - y(ii,jj))^2;
%             s(i,2) = s(i,1) + (ii - jj)* sign(x2(ii,jj) - y(ii,jj))*(x2(ii,jj) - y(ii,jj))^2;
%             s(i,3) = s(i,1) + (ii - jj)* sign(x3(ii,jj) - y(ii,jj))*(x3(ii,jj) - y(ii,jj))^2;
%         end
%     end


% figure
% plot (g);
% legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');
% set(gca,'XTick',1:1:4);
% set(gca, 'XTickLabel', {'09Q1','09Q2','09Q3','09Q4'});
% title('average Rating Default Probability Forecast: 2009Q1 - 2009Q4');
%
%
% figure
% plot (s);
% legend('model','benchmark1-avg','benchmark2-previous','Location','NorthWest');
% set(gca,'XTick',1:1:4);
% set(gca, 'XTickLabel', {'09Q1','09Q2','09Q3','09Q4'});
% title('D3 Compariaon: Forecast and Real: 2009Q1 - 2009Q4');

% % visual check est. TM-----------------------------------------------------
%
% % first, we look at criteria D3
% g1 (:,1)= Pfcst(i,1,1);  g1(:,2) = naive1(:,1,1); g1(:,3) = naive2(:,1,1); g1(:,4)= stdProbQ (1:length(PHat),1,1);% rank 1.5 - 1.5
% g2 (:,1)= PHat(8,:,3,3); g2(:,2) = naive1(:,3,3); g2(:,3) = naive2(:,3,3); g2(:,4)= stdProbQ (1:length(PHat),3,3);% rank 2.5 - 2.5
% g3 (:,1)= PHat(8,:,6,6); g3(:,2) = naive1(:,6,6); g3(:,3) = naive2(:,6,6); g3(:,4)= stdProbQ (1:length(PHat),6,6);% rank 3 - 3
% g4 (:,1)= PHat(8,:,5,8); g4(:,2) = naive1(:,5,8); g4(:,3) = naive2(:,5,8); g4(:,4)= stdProbQ (1:length(PHat),5,8);% rank 3.5 -D
% g5 (:,1)= PHat(8,:,6,8); g5(:,2) = naive1(:,6,8); g5(:,3) = naive2(:,6,8); g5(:,4)= stdProbQ (1:length(PHat),6,8);% rank 4 -D
% g6 (:,1)= PHat(8,:,7,8); g6(:,2) = naive1(:,7,8); g6(:,3) = naive2(:,7,8); g6(:,4)= stdProbQ (1:length(PHat),7,8);% rank 4.5 -D
%
% % then, we look at criteria SVD
%
% g6 (:,1)= PHat(5,:,7,8); g6(:,2) = naive1(:,7,8); g6(:,3) = naive2(:,7,8); g6(:,4)= stdProbQ (1:length(PHat),7,8);% rank 4.5 -D
%
% figure
% plot (g6); %([g1,gg1,ggg1]);%,'-*'); plot (gg1,'-rs');
% set(gca,'XTick',1:4:length(g1));
% labels = quaterlabels(1998, length(g1));
% labels = labels(1:4:length(g1));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');




% % W = zeros (outSampleLength, 3); % first column is for metric BFS, the second colum is for metric D2
% PHat_BFS = zeros (inSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% % PHat_D2 = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% % PHat_SVD = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% R = RATE_LIST_LENGTH-1;
% C = RATE_LIST_LENGTH;
% yrEst= INSAMPLE_YR - 1998 + 1 + floor(INSAMPLE_QR / 4);
%
% % for i = 1: inSampleLength %outSampleLength
%
% P = stdProbQ;
% X(:,:) = bin(inSampleLength,:,:);
% wantedZ = Z;
% t = inSampleLength;
%
% % BFS distance
% obj1=@(w)difBFS_fixedW(w,t,N,P,X,wantedZ,R,C);
% [foundW,fval,exitflag,output] = fminbnd(obj1,0,0.5);
% figure;
% fplot(obj1,[-0.5 0.5],'*')
% title('Mininization Object with BFS distance',...
%   'FontWeight','bold')
% PHat_BFS = calculatePHat_fixedW(foundW,t,X,wantedZ,R,C);
%
%
% % W(i,1) = foundW;
% % sheetlist=strcat('',year_list_text(yrEst,:),'.',quarter_list_text(qrEst,:),'');
% % xlswrite(filename, PHat, sheetlist);
% % figure;
% % subplot(outSampleLength,1,i); fplot(obj1,[-0.5 0.5],'*')
% % title('Mininization Object with BFS distance',...
% %   'FontWeight','bold')
%
%
% % D2 distance
% obj2=@(w)difD2(w,P,X,wantedZ,R,C);
% [foundW,fval,exitflag,output] = fminbnd(obj2,0,0.5);
% W(i,2) = foundW;
% PHat_D2(i,:,:) = calculatePHat(foundW,X,wantedZ,R,C);
%
%
% % SVD
% obj3=@(w)difSVD_fixedW(w,t,P,X,wantedZ,R,C);
% [foundW,fval,exitflag,output] = fminbnd(obj3,0,0.5);
% % W(i,3) = foundW;
% PHat_SVD = calculatePHat_fixedW(foundW,t,X,wantedZ,R,C);
% % xx(:,:) = PHat_SVD(i,:,:);
% % sheetlist=strcat('',year_list_text(yrEst,:),'.',quarter_list_text(qrEst,:),'');
% % xlswrite(filename, xx, sheetlist);
% % clear xx;
% % end




% MAE ---------------------------------------------------------------------
% mae_perf = zeros(5,outSampleLength);
% mae_svd_perf = zeros(5,outSampleLength);
%
% for i = 1: outSampleLength
%     P(:,:) = stdProbQ(inSampleLength+i,:,:);
%     previous (:,:) = stdProbQ(inSampleLength+i-1,:,:);
%     average (:,:) = avgRollingTM(inSampleLength+i-1,:,:);
%     model_BFS (:,:) = PHat_BFS(i,:,:);
%     model_D2 (:,:) = PHat_D2(i,:,:);
%     model_SVD (:,:) = PHat_SVD(i,:,:);
%
%
%     mae_perf(1,i) = mean2 (abs (P - previous));
%     mae_perf(2,i) = mae (P - average);
%     mae_perf(3,i) = mae (P - model_BFS);
%     mae_perf(4,i) = mae (P - model_D2);
%     mae_perf(5,i) = mae (P - model_SVD);
%
%     mae_svd_perf(1,i) =  abs (mean(svd(P))- mean(svd(previous)));
%     mae_svd_perf(2,i) =  abs (mean(svd(P))- mean(svd(average)));
%     mae_svd_perf(3,i) =  abs (mean(svd(P))- mean(svd(model_BFS)));
%     mae_svd_perf(4,i) =  abs (mean(svd(P))- mean(svd(model_D2)));
%     mae_svd_perf(5,i) =  abs (mean(svd(P))- mean(svd(model_SVD)));
%
%
% end

% naive1=stdProbAvgQ;
% naive2=stdProbQ(inSampleLength+i,:,:);
% % n(:,:)=N(yrp,freqjp,2:10,13);
% % nn=repmat(n,1,10);
% % nt(1,1,:,:)=nn;
% % Preal=P(:,:,2:10,2:11);
% e_mod=(nt.*(Preal-PZ).^2)./(PZ.*(I-PZ));
% e_avg=(nt.*(Preal-naive1).^2)./(naive1.*(I-naive1));
% e_pre=(nt.*(Preal-naive2).^2)./(naive2.*(I-naive2));
% perf_mod = mae(e_mod);
% perf_naive1 = mae(e_avg);
% perf_naive2 = mae(e_pre);
% %
%
% Parameters={'foundW','Zhat', 'perf_mod', 'perf_naive1','perf_naive2'; foundW Zhat perf_mod perf_naive1 perf_naive2}';
% xlswrite(filename, Parameters,'Parameters');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Naive1 MAE of L1,L2 NSD,D1,D2   %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% L1_n1=Preal-naive1;
% L2_n1=(Preal-naive1).^2;
% NSD_n1=((Preal-naive1).^2)./(naive1);
% % D1_n1=
% % D2_n1=
%
% perf_naive1e_L1=mae(L1_n1);
% perf_naive1e_L2=sqrt(mae(L2_n1).*90)/90;
% perf_naive1e_NSD=mae(NSD_n1);
% % perf_naive1e_D1=
% % perf_naive1e_D2=
%
% disp('fini');




% for i = 1 : outSampleLength
%
% % SVDhat = SVD (1:inSampleLength + i -1, :);
% % SVDhat(inSampleLength + i, :) = fcstSVD (i, :);
% % ZhatVector = - zscore (SVDhat);
% % Zhat(i, :) = ZhatVector(size(ZhatVector,1),:);
% % save in 'optimization_07q4_jun2_3.xls';
%
% Zhat(i,:) = - ( fcstInvPD(i,:)- mean(inversePD(1:inSampleLength+i-1,:)) )/std(inversePD(1:inSampleLength+i-1,:));
% % save in 'optimization_07q4_jun2_5.xls';
%
% end

